---
title: "CE Index Analysis File"
output: html_document
date: "2025-02-07"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
rm(list = ls())
cat("\014")


library(tidyverse)
library(vroom)
library(ggrepel)
library(ggthemes)
library(forcats)
library(fixest)
library(zoo)

options(scipen = 99)
```

*Load Data* 

```{r}
test <- vroom("final.data.csv")


test <- test %>%
  distinct(state, year, .keep_all = T) %>%
  filter(!state %in% c("US", "DC")) %>%
  mutate(state_id = as.numeric(factor(state))) %>%
  filter(year > 1999, year < 2023)


pol <- vroom("data/final.partisan.data.csv")
fips <- vroom("data/FIPS_CODES.csv")

pol <- pol %>%
  mutate(statename = tolower(state)) %>%
  dplyr::select(-...1) %>%
  left_join(fips) %>%
  dplyr::select(stateabbr, year, two.party.dem.senate.share, 
                two.party.dem.house.share) %>%
  rename(state = stateabbr) %>%
  distinct(state, year, .keep_all = T)


test <- left_join(test, pol, by = c("state", "year"))

test <- test %>%
  distinct(state, year, .keep_all = T) %>%
  filter(year >= 1995 & year < 2024)

df <- vroom("data/states.data.csv") %>%
  dplyr::select(abb, year, bea_region, climate) %>%
  rename(state = abb)

test <- left_join(test, df, by = c("state", "year"))

test <- test %>%
  group_by(state) %>%
  mutate(bea_region = ifelse(is.na(bea_region), unique(bea_region[!is.na(bea_region)]), bea_region)) %>%
  ungroup()

test <- test %>%
  distinct(state, year, .keep_all = T)


test <- test %>%
  dplyr::select(state, state_id, year, Latent_Factor, climate, two.party.dem.senate.share,
                two.party.dem.house.share, tot.energy.bn.btu.cons,
                perc.renew.energy.prod.bn.btu, perc.renew.energy.cons,
                energy.cons.per.2017.usd.gdp, local.energy.exps.pc,
                local.energy.price, mt.co2.per.mn.2017.gdp, ff.gdp.perc,
                gdp.mn.2017.usd, th.mt.mat.cons.per.mn.2017.gdp)
```



```{r}
gdp.data <- vroom("data/industry.gdp.data.csv") %>%
  dplyr::select(state, year, pop.in.th, manuf.gdp.perc)

test <- left_join(test, gdp.data, by = c("state", "year"))

```

```{r}
rm(fips, gdp.data, pol)
```


*Descriptives Clean Economy Index*

```{r}
data_2000 <- test %>%
  filter(year == 2000) %>%
  mutate(mean = mean(Latent_Factor),
         sd = sd(Latent_Factor)) %>%
  dplyr::select(mean, sd)

data_2022 <- test %>%
  filter(year == 2022) %>%
  mutate(mean = mean(Latent_Factor),
         sd = sd(Latent_Factor)) %>%
  dplyr::select(mean, sd)
```


*Comparison with SDG progress assessments*

```{r}
sdr <- read.csv("data/SDR_2021.csv")

sdg_2018 <- read.csv("data/sdg_2018.csv")
```

```{r}
data_2021 <- test %>%
  filter(year == 2021) %>%
  left_join(sdr) %>%
  dplyr::select(state, year, Latent_Factor, sdg7_score, sdg13_score, sdg12_score) %>%
  mutate(year = 2021)

data_2018 <- test %>%
  filter(year == 2018) %>%
  left_join(sdg_2018, by = c("state", "year")) %>%
  dplyr::select(state, year, Latent_Factor, SDG07, SDG13, SDG12) %>%
  rename(sdg7_score = SDG07, sdg13_score = SDG13, sdg12_score = SDG12)

cor_data <- rbind(data_2018, data_2021)

```

*Pooled 2021 and 2018* 
```{r}

cor(cor_data$Latent_Factor, cor_data$sdg7_score, method = "pearson")

cor(cor_data$Latent_Factor, cor_data$sdg12_score, method = "pearson")

cor(cor_data$Latent_Factor, cor_data$sdg13_score, method = "pearson")


cor(cor_data$Latent_Factor, cor_data$sdg7_score, method = "spearman")

cor(cor_data$Latent_Factor, cor_data$sdg12_score, method = "spearman")

cor(cor_data$Latent_Factor, cor_data$sdg13_score, method = "spearman")

```

```{r}
#Only 2021

cor(data_2021$Latent_Factor, data_2021$sdg7_score, method = "pearson")

cor(data_2021$Latent_Factor, data_2021$sdg12_score, method = "pearson")

cor(data_2021$Latent_Factor, data_2021$sdg13_score, method = "pearson")


cor(data_2021$Latent_Factor, data_2021$sdg7_score, method = "spearman")

cor(data_2021$Latent_Factor, data_2021$sdg12_score, method = "spearman")

cor(data_2021$Latent_Factor, data_2021$sdg13_score, method = "spearman")
```

```{r}
# Only 2018 

cor(data_2018$Latent_Factor, data_2018$sdg7_score, method = "pearson")

cor(data_2018$Latent_Factor, data_2018$sdg12_score, method = "pearson")

cor(data_2018$Latent_Factor, data_2018$sdg13_score, method = "pearson")


cor(data_2018$Latent_Factor, data_2018$sdg13_score, method = "spearman")

cor(data_2018$Latent_Factor, data_2018$sdg7_score, method = "spearman")

cor(data_2018$Latent_Factor, data_2018$sdg12_score, method = "spearman")
```


*Visualization of scores by state*
```{r}

test %>%
  filter(year == 2022) %>%
  arrange(desc(Latent_Factor)) %>%
  ggplot(aes(x = reorder(state, Latent_Factor), y = Latent_Factor)) +
  geom_point(size = 0.8) +
  coord_flip() +
  theme_minimal() +
  xlab("") +
  ylab("") +
  scale_y_continuous(expand = expansion(mult = c(0.1, 0.1))) + 
  theme(axis.text.y = element_text(size = 6, hjust = 1, vjust = 0.5))


ggsave("cei.final.png", 
       width = 8, height = 8/1.618, units = "in")
```

*Visualization of scores over time*

```{r}
test %>%
  filter(year < 2023 & year > 1999) %>%
 ggplot(aes(x = year, y = Latent_Factor, group = state)) +
  geom_line(alpha = 0.3, size = 0.7) + 
  theme_clean() + 
  xlab("") + ylab("Green Transition Index")
```


```{r}

highlight_states <- c("CA", "WY", "TX")

test %>%
  filter(year > 1996, year < 2023) %>%

  ggplot(aes(x = year, y = Latent_Factor, group = state)) +
  geom_line(
    data = . %>% filter(!state %in% highlight_states),
    color = "grey80", alpha = 0.5, size = 0.7
  ) +
  geom_line(
    data = . %>% filter(state %in% highlight_states),
    aes(color = state),
    size = 1.2
  ) +
  scale_color_manual(
    values = c(
      CA = "#1f78b4",  
      WY = "#e31a1c",  
      TX = "#33a02c"   
    ),
    name = "State"
  ) +

  theme_clean() +
  xlab(NULL) +
  ylab("Green Transition Index") +
  theme(
    legend.position = "bottom"
  )

ggsave("lineplot.final.png", 
       width = 8, height = 8/1.618, units = "in")
```


**Scatterplot**

```{r}
test %>%
  filter(year == 2019) %>%
  ggplot(aes(x = climate, y = Latent_Factor)) +
  geom_point(aes(color = two.party.dem.house.share)) + 
  geom_text_repel(aes(label = state, color = two.party.dem.house.share)) + 
  theme_clean() + 
  scale_color_gradient2(
    low = "red",
    mid = "grey",
    high = "blue",
    midpoint = 50 
  ) +
  geom_vline(aes(xintercept = median(climate, na.rm = TRUE)), linetype = "dashed", color = "black") + 
  geom_hline(aes(yintercept = median(Latent_Factor, na.rm = TRUE)), linetype = "dashed", color = "black") + 
  guides(color = "none") + 
  labs(
    title = "",
    x = "Climate Policy Stringency",
    y = "Green Transition Index",
    caption = ""
  )


ggsave("scatter.color.png", 
       width = 8, height = 8/1.618, units = "in")
```

**Map** 

```{r}
# Rolling 3-year average

test <- test %>%
  group_by(state) %>% 
  arrange(year) %>%    
  mutate(
    rolling_avg = rollmean(Latent_Factor, k = 3, fill = NA, align = "center") 
  )

data <- test %>%
  dplyr::select(state, year, Latent_Factor, rolling_avg) %>%
  filter(year > 1999 & year < 2024) %>%
  filter(!state %in% c("DC", "US"))
```


```{r}
library(maps)
library(usmap)

data_2020 <- data %>%
  filter(year == 2021) %>%
  mutate(state = tolower(state)) 


plot_usmap(data = data_2020, values = "rolling_avg", regions = "states") +
  scale_fill_gradient2(
    low = "red",       
    mid = "white",     
    high = "green",    
    midpoint = 0,      
    na.value = "grey"  
  ) +
  labs(
    title = "Rolling 3-Year Average Across US States (2021)",
    subtitle = "",
    fill = "Rolling Avg"
  ) +
  theme(
    panel.background = element_rect(color = "black", fill = "white"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 7),  
    legend.key.size = unit(0.2, "cm"),    
    legend.key.width = unit(0.9, "cm"),
    legend.direction = "horizontal"
  )

width <- 8 

ggsave("Map.2021.png", width = width, height = width/1.618, units = "in",
       dpi = 600)
```


*Combined map* 
```{r}

data_2000 <- data %>%
  filter(year == 2001) %>%
  mutate(state = tolower(state)) 


combined_range <- range(c(data_2000$rolling_avg, data_2020$rolling_avg), na.rm = TRUE)
print(combined_range)

```

```{r}
library(gridExtra)

# Map for 2000
map_2000 <- plot_usmap(data = data_2000, values = "rolling_avg", regions = "states") +
  scale_fill_gradient2(
    low = "red",
    mid = "white",
    high = "green",
    midpoint = 0,
    na.value = "grey",
    limits = combined_range  
  ) +
  labs(
    title = "2000",
    fill = "Rolling Avg"
  ) +
  theme(
    panel.background = element_rect(color = "black", fill = "white"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 7),
    legend.key.size = unit(0.2, "cm"),
    legend.key.width = unit(0.9, "cm"),
    legend.direction = "horizontal"
  )

# Map for 2020
map_2020 <- plot_usmap(data = data_2020, values = "rolling_avg", regions = "states") +
  scale_fill_gradient2(
    low = "red",
    mid = "white",
    high = "green",
    midpoint = 0,
    na.value = "grey",
    limits = combined_range  
  ) +
  labs(
    title = "2020",
    fill = "Rolling Avg"
  ) +
  theme(
    panel.background = element_rect(color = "black", fill = "white"),
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 7),
    legend.key.size = unit(0.2, "cm"),
    legend.key.width = unit(0.9, "cm"),
    legend.direction = "horizontal"
  )

library(patchwork)


map_2000 <- map_2000 + theme(legend.position = "none")
map_2020 <- map_2020 + theme(legend.position = "none")


combined_plot <- (map_2000 | map_2020) +
  plot_annotation(
    title = "",
    theme = theme(plot.title = element_text(size = 14, face = "bold"))
  ) &
  theme(
    legend.position = "bottom"
  )


combined_plot <- combined_plot & 
  guides(fill = guide_colorbar(title = "Rolling Avg"))


print(combined_plot)


ggsave("Map.2000.2023.png", width = width, height = width/1.618, units = "in",
       dpi = 600)

```


**Analysis** 

```{r}
library(sandwich)
library(lmtest)
library(plm)
library(estimatr)
```


*TWFE Model with Climate Policy Stringency  from Berquist and Warshaw (2023)*


*Lagged variables of climate policy*

```{r}
df <- test %>%
  group_by(state) %>%
  arrange(state, year) %>%
  mutate(climate_lag1 = dplyr::lag(climate, n = 1)) %>%
  mutate(climate_lag2 = dplyr::lag(climate, n = 2)) %>%
  mutate(DV_lag1 = dplyr::lag(Latent_Factor, n = 1)) %>%
  mutate(DV_lag2 = dplyr::lag(Latent_Factor, n = 2)) %>%
  mutate(DV_lag3 = dplyr::lag(Latent_Factor, n = 3)) %>%
  ungroup()

```


*Dynamic TWFE w 1lag* 

*One lag*
```{r}
reg <- feols(Latent_Factor ~ climate + DV_lag1 | state + year, data = df, cluster = ~state)

summary(reg)

BIC(reg)
```

*1-year lag*
```{r}
reg <- feols(Latent_Factor ~ climate_lag1 + DV_lag1 | state + year, data = df, cluster = ~state)

summary(reg)

```

*2-year lag*
```{r}
reg <- feols(Latent_Factor ~ climate_lag2 + DV_lag1 | state + year, data = df, cluster = ~state)

summary(reg)

```


*Simple TWFE*

```{r}
reg <- feols(Latent_Factor ~ climate | state + year, data = df, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Latent_Factor ~ climate_lag2 | state + year, data = df, cluster = ~state)

summary(reg)
```




*Alternative measure of Climate Policy - RPS schemes*

*RPS as alternative Climate Policy Measure*

```{r}
rps <- vroom("data/RPS.csv")
rps <- rps %>%
  dplyr::select(1:25) %>%
  pivot_longer(cols = 3:25, names_to = "year", values_to = "RPS_perc") %>%
  dplyr::select(State, year, RPS_perc) %>%
  rename(state = State) %>%
  mutate(year = as.numeric(year))

rps$RPS_perc <- as.numeric(gsub("%", "", rps$RPS_perc))

df <- left_join(test, rps, by = c("state", "year"))


df$RPS_perc[is.na(df$RPS_perc)] <- 0

df <- df %>%
  mutate(any_rps = if_else(RPS_perc > 0, 1,0))
```

*Lagged RPS*

```{r}
df <- df %>%
  group_by(state) %>%
  arrange(state, year) %>%
  mutate(rps_lag2 = dplyr::lag(any_rps, n = 2)) %>%
  mutate(rps_perc_lag2 = dplyr::lag(RPS_perc, n = 2)) %>%
   mutate(DV_lag1 = dplyr::lag(Latent_Factor, n = 1)) %>%
  ungroup()
```

*Simple TWFE* 
```{r}
reg <- feols(Latent_Factor ~ rps_lag2 | state + year, data = df, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Latent_Factor ~ rps_lag2 + DV_lag1| state + year, data = df, cluster = ~state)

summary(reg)

```


**Callaway and Sant'anna -- RPS dummies** 

```{r}
treat <- df %>%
  group_by(state) %>%
  mutate(first_treated = if_else(rps_lag2 == 1, year, NA_integer_)) %>%
  summarize(first_treated = min(first_treated, na.rm = TRUE)) %>%
  mutate(first_treated = if_else(is.infinite(first_treated), 0, first_treated))

test <- left_join(data, treat)

library(did)

did_results <- att_gt(
  yname = "Latent_Factor",         
  tname = "year",       
  idname = "state_id",     
  gname = "first_treated", 
  xformla = ~ 1,        
  data = test,
  panel = FALSE,
  clustervars = "state_id"       

)
```


```{r}
agg.simple <- aggte(did_results, type = "simple")

mw.dyn <- aggte(did_results, type = "dynamic")

agg.gs <- aggte(did_results, type = "group")

```


```{r}

results <- list(
  Group_Time = agg.simple,
  Dynamic = mw.dyn,
  Group = agg.gs
)

pe.vis <- do.call(rbind, lapply(names(results), function(name) {
  obj <- results[[name]]
  data.frame(
    Treatment = name,
    ATT = obj$overall.att,
    SE = obj$overall.se,
    Lower_CI = obj$overall.att - (1.96 * obj$overall.se),
    Upper_CI = obj$overall.att + (1.96 * obj$overall.se)
  )
}))

ggplot(pe.vis, aes(x = Treatment, y = ATT, ymin = Lower_CI, ymax = Upper_CI)) +
  geom_point(linewidth = 2, color = "black") +
  geom_linerange(aes(ymin = Lower_CI, ymax = Upper_CI), color = "black", size = 0.5, alpha = 0.7) +
  labs(
    title = "Treatment Effect of RPS by Aggregation Method",
    x = "",
    y = ""
  ) +
  theme_classic() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  coord_flip()

ggsave("Treatment_effect_Callaway_SantAnna_Estimator-RPS.png", 
       width = 8, height = 8/1.618, units = "in")

```



**Partisanship**

*Exclude Nebraska*
```{r}
df <- df %>%
  filter(!state == "NE")
```

*Create Republican control* 
```{r}
df <- df %>%
  mutate(rep_control = if_else(two.party.dem.house.share < 50 & two.party.dem.senate.share < 50, 1, 0), 
         rep_part = if_else(two.party.dem.house.share < 50 | two.party.dem.senate.share < 50, 1, 0))

```


*Create level of Republican control as ordinal variable* 

```{r}

govs <- vroom("data/governors.csv")

data <- left_join(df, govs, by = c("state", "year"))

# 
data <- data %>%
  mutate(rep_gov = if_else(gov_party == "Rep", 1, 0))


data <- data %>%
  mutate(
    unified_dem = as.integer(
      two.party.dem.house.share > 50 & 
      two.party.dem.senate.share > 50 & 
      gov_party == "Dem"
    ),
    unified_rep = as.integer(
      two.party.dem.house.share < 50 & 
      two.party.dem.senate.share < 50 & 
      gov_party == "Rep"
    ),
    divided = as.integer(
      !(unified_dem == 1 | unified_rep == 1)
    )
  )

```


```{r}
data <- data %>%
  group_by(state) %>%
  arrange(state, year) %>%
  mutate(rep_gov_lag2 = dplyr::lag(rep_gov, n = 2)) %>%
  mutate(rep_control_lag2 = dplyr::lag(rep_control, n = 2)) %>%
  mutate(unified_rep_lag1 = dplyr::lag(unified_rep, n = 1)) %>%
  mutate(unified_rep_lag2 = dplyr::lag(unified_rep, n = 2)) %>%
  mutate(unified_dem_lag1 = dplyr::lag(unified_dem, n = 1)) %>%
  mutate(unified_dem_lag2 = dplyr::lag(unified_dem, n = 2)) %>%
    mutate(climate_lag2 = dplyr::lag(climate, n = 2)) %>%
   mutate(DV_lag1 = dplyr::lag(Latent_Factor, n = 1)) %>%
  ungroup()
```



*Unified Rep Control*

*Dynamic*
```{r}
reg <- feols(Latent_Factor ~ unified_rep + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Latent_Factor ~ unified_rep_lag1 + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```

```{r}
reg <- feols(Latent_Factor ~ unified_rep_lag2 + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```


*Simple*

```{r}
reg <- feols(Latent_Factor ~ unified_rep | state + year, data = data, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Latent_Factor ~ unified_rep_lag2  | state + year, data = data, cluster = ~state)

summary(reg)
```




*Unified Dem Control*

```{r}
reg <- feols(Latent_Factor ~ unified_dem + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```

*Lagged*

```{r}
reg <- feols(Latent_Factor ~ unified_dem_lag1 + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Latent_Factor ~ unified_dem_lag2 + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```

*Simple TWFE*

```{r}
reg <- feols(Latent_Factor ~ unified_dem_lag2 | state + year, data = data, cluster = ~state)

summary(reg)
```

```{r}
reg <- feols(Latent_Factor ~ unified_dem | state + year, data = data, cluster = ~state)

summary(reg)
```



**Callaway and Sant'anna -- Robustness test for Partisanship dummies** 


```{r}
treat <- data %>%
  group_by(state) %>%
  mutate(first_treated = if_else(unified_rep == 1, year, NA_integer_)) %>%
  summarize(first_treated = min(first_treated, na.rm = TRUE)) %>%
  mutate(first_treated = if_else(is.infinite(first_treated), 0, first_treated))

test <- left_join(data, treat)

library(did)

did_results <- att_gt(
  yname = "Latent_Factor",         
  tname = "year",       
  idname = "state_id",     
  gname = "first_treated", 
  xformla = ~ 1,        
  data = test,
  panel = FALSE,
  

)
```

```{r}
agg.simple <- aggte(did_results, type = "simple")

mw.dyn <- aggte(did_results, type = "dynamic")

agg.gs <- aggte(did_results, type = "group")

```


```{r}
results <- list(
  Group_Time = agg.simple,
  Dynamic = mw.dyn,
  Group = agg.gs
)

pe.vis <- do.call(rbind, lapply(names(results), function(name) {
  obj <- results[[name]]
  data.frame(
    Treatment = name,
    ATT = obj$overall.att,
    SE = obj$overall.se,
    Lower_CI = obj$overall.att - (1.96 * obj$overall.se),
    Upper_CI = obj$overall.att + (1.96 * obj$overall.se)
  )
}))

ggplot(pe.vis, aes(x = Treatment, y = ATT, ymin = Lower_CI, ymax = Upper_CI)) +
  geom_point(linewidth = 2, color = "black") +
  geom_linerange(aes(ymin = Lower_CI, ymax = Upper_CI), color = "black", size = 0.5, alpha = 0.7) +
  labs(
    title = "Treatment Effect of Republican State Control by Aggregation Method",
    x = "",
    y = ""
  ) +
  theme_classic() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  coord_flip()

# Save the plot
ggsave("Treatment_effect_Callaway_SantAnna_Estimator-1.png", 
       width = 8, height = 8/1.618, units = "in")

```


```{r}

test <- test %>%
  group_by(state_id) %>%
  mutate(state_id_treatment = as.numeric(factor(paste(state_id, rep_control, sep = "_"))))



did_results <- att_gt(
  yname = "Latent_Factor",         
  tname = "year",                  
  idname = "state_id_treatment",   
  gname = "first_treated",         
  xformla = ~ 1,                   
  data = test,
  panel = FALSE
)


agg.simple <- aggte(did_results, type = "simple")

mw.dyn <- aggte(did_results, type = "dynamic")

agg.gs <- aggte(did_results, type = "group")

results <- list(
  Group_Time = agg.simple,
  Dynamic = mw.dyn,
  Group = agg.gs
)


pe.vis <- do.call(rbind, lapply(names(results), function(name) {
  obj <- results[[name]]
  data.frame(
    Treatment = name,
    ATT = obj$overall.att,
    SE = obj$overall.se,
    Lower_CI = obj$overall.att - (1.96 * obj$overall.se),
    Upper_CI = obj$overall.att + (1.96 * obj$overall.se)
  )
}))


ggplot(pe.vis, aes(x = Treatment, y = ATT, ymin = Lower_CI, ymax = Upper_CI)) +
  geom_point(linewidth = 2, color = "black") +
  geom_linerange(aes(ymin = Lower_CI, ymax = Upper_CI), color = "black", size = 0.5, alpha = 0.7) +
  labs(
    title = "Treatment Effect of Republican State Control by Aggregation Method",
    x = "",
    y = ""
  ) +
  theme_classic() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  coord_flip()


ggsave("Treatment_effect_Callaway_SantAnna_Estimator-2.png", 
       width = 8, height = 8/1.618, units = "in")
```



*Alternative measures of Partisanship*

```{r}
data <- data %>%
  group_by(state) %>%
  arrange(state, year) %>%
  mutate(rep_gov_lag2 = dplyr::lag(rep_gov, n = 2)) %>%
  mutate(rep_control_lag2 = dplyr::lag(rep_control, n = 2)) %>%
  mutate(rep_part_lag2 = dplyr::lag(rep_part, n = 2)) %>%
  mutate(unified_rep_lag2 = dplyr::lag(unified_rep, n = 2)) %>%
  mutate(unified_dem_lag2 = dplyr::lag(unified_dem, n = 2)) %>%
    mutate(climate_lag2 = dplyr::lag(climate, n = 2)) %>%
   mutate(DV_lag1 = dplyr::lag(Latent_Factor, n = 1)) %>%
  ungroup()
```


**Republican Control of lower-house** 

```{r}
reg <- feols(Latent_Factor ~ rep_control + DV_lag1| state + year, data = data, cluster = ~state)

summary(reg)

```

*Any Republican majority in any house*

```{r}
reg <- feols(Latent_Factor ~ rep_part +DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```


*Rep Governors*

```{r}
reg <- feols(Latent_Factor ~ rep_gov + DV_lag1 | state + year, data = data, cluster = ~state)

summary(reg)
```



*Robustness Test*

**Core index w fewer inputs** 

```{r}
core <- read.csv("core.index.csv")

core_climate <- left_join(test, core, by = c("state", "year")) %>%
  group_by(state) %>%
  arrange(state, year) %>%
  mutate(climate_lag1 = dplyr::lag(climate, n = 1)) %>%
  mutate(DV_lag1 = dplyr::lag(Core_index, n = 1))


sd <- data %>%
  dplyr::select(state, year, Latent_Factor, Core_index) %>%
  filter(year < 2023, year > 1999)

cor(sd$Latent_Factor, sd$Core_index, method = "pearson")
```

*Dynamic* 

```{r}
reg <- feols(Core_index ~ climate + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)

summary(reg)
```

*Dynamic w lag * 

```{r}
reg <- feols(Core_index ~ climate_lag1 + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)

summary(reg)
```



*Simple*

```{r}
reg <- feols(Core_index ~ climate | state + year, data = core_climate, 
             cluster = ~state)

summary(reg)
```

*Simple w lag*

```{r}
reg <- feols(Core_index ~ climate_lag1 | state + year, data = core_climate, 
             cluster = ~state)

summary(reg)
```



*Dynamic RPS*

```{r}

core_climate <- left_join(test, core, by = c("state", "year")) %>%
  group_by(state) %>%
  arrange(state, year) %>%
  mutate(climate_lag1 = dplyr::lag(climate, n = 1)) %>%
  mutate(DV_lag1 = dplyr::lag(Core_index, n = 1))


rps_core <- df %>%
  dplyr::select(state, year, any_rps)

rps_core <- left_join(core_climate, rps_core, by = c("state", "year"))

rps_core <- rps_core %>%
  mutate(rps_lag1 = dplyr::lag(any_rps))


reg <- feols(Core_index ~ any_rps + DV_lag1 | state + year, data = rps_core, 
             cluster = ~state)

summary(reg)
```


*Simple*

```{r}
reg <- feols(Core_index ~ rps_lag1 + DV_lag1 | state + year, data = rps_core, 
             cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Core_index ~ any_rps| state + year, data = rps_core, 
             cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Core_index ~ rps_lag1 | state + year, data = rps_core, 
             cluster = ~state)

summary(reg)
```




Republican control

*Dynamic*

```{r}
robust <- left_join(data, core, by = c("state", "year"))

reg <- feols(Core_index ~ unified_rep + DV_lag1 | state + year, data = robust, cluster = ~state)

summary(reg)
```

```{r}
reg <- feols(Core_index ~ unified_rep_lag1 + DV_lag1| state + year, data = robust, cluster = ~state)

summary(reg)
```




*Simple*
```{r}
reg <- feols(Core_index ~ unified_rep | state + year, data = robust, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Core_index ~ unified_rep_lag1 | state + year, data = robust, cluster = ~state)

summary(reg)
```



*Unified Dem Control*

*dynamic*
```{r}
reg <- feols(Core_index ~ unified_dem + DV_lag1 | state + year, data = robust, cluster = ~state)

summary(reg)
```


```{r}
reg <- feols(Core_index ~ unified_dem_lag1 + DV_lag1 | state + year, data = robust, cluster = ~state)

summary(reg)
```


*Simple*
```{r}
reg <- feols(Core_index ~ unified_dem | state + year, data = robust, cluster = ~state)

summary(reg)

```


```{r}
reg <- feols(Core_index ~ unified_dem_lag1 | state + year, data = robust, cluster = ~state)

summary(reg)

```





**Robustness to running model on key individual indicators**

*Climate Policy*

```{r}
core_climate <- core_climate %>%
  mutate(DV_lag1 = dplyr::lag(perc.renew.energy.prod.bn.btu, n = 1))

reg <- feols(perc.renew.energy.prod.bn.btu ~ climate + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)

summary(reg)
```

```{r}

core_climate <- core_climate %>%
  mutate(DV_lag1 = dplyr::lag(perc.renew.energy.cons, n = 1))

reg <- feols(perc.renew.energy.cons ~ climate + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)
summary(reg)


```

```{r}
core_climate <- core_climate %>%
  mutate(DV_lag1 = dplyr::lag(energy.cons.per.2017.usd.gdp, n = 1))

reg <- feols(energy.cons.per.2017.usd.gdp ~ climate + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)
summary(reg)
```

```{r}
core_climate <- core_climate %>%
  mutate(DV_lag1 = dplyr::lag(mt.co2.per.mn.2017.gdp, n = 1))

reg <- feols(mt.co2.per.mn.2017.gdp ~ climate + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)
summary(reg)

```

```{r}
core_climate <- core_climate %>%
  mutate(DV_lag1 = dplyr::lag(th.mt.mat.cons.per.mn.2017.gdp, n = 1))

reg <- feols(th.mt.mat.cons.per.mn.2017.gdp ~ climate + DV_lag1 | state + year, data = core_climate, 
             cluster = ~state)
summary(reg)
```

```{r}
core_climate <- core_climate %>%
  mutate(DV_lag1 = dplyr::lag(ff.gdp.perc, n = 1))

reg <- feols(ff.gdp.perc ~ climate + DV_lag1| state + year, data = core_climate, 
             cluster = ~state)
summary(reg)
```



*Republicans*

```{r}

```

```{r}
robust <- robust %>%
  mutate(DV_lag1 = dplyr::lag(perc.renew.energy.prod.bn.btu, n = 1))

reg <- feols(perc.renew.energy.prod.bn.btu ~ unified_rep + DV_lag1| state + year, data = robust, 
             cluster = ~state)

summary(reg)
```

```{r}
robust <- robust %>%
  mutate(DV_lag1 = dplyr::lag(perc.renew.energy.cons, n = 1))

reg <- feols(perc.renew.energy.cons ~ unified_rep + DV_lag1| state + year, data = robust, 
             cluster = ~state)
summary(reg)


```

```{r}
robust <- robust %>%
  mutate(DV_lag1 = dplyr::lag(energy.cons.per.2017.usd.gdp, n = 1))

reg <- feols(energy.cons.per.2017.usd.gdp ~ unified_rep + DV_lag1| state + year, data = robust, 
             cluster = ~state)
summary(reg)
```

```{r}
robust <- robust %>%
  mutate(DV_lag1 = dplyr::lag(mt.co2.per.mn.2017.gdp, n = 1))

reg <- feols(mt.co2.per.mn.2017.gdp ~ unified_rep + DV_lag1| state + year, data = robust, 
             cluster = ~state)
summary(reg)

```

```{r}
robust <- robust %>%
  mutate(DV_lag1 = dplyr::lag(th.mt.mat.cons.per.mn.2017.gdp, n = 1))

reg <- feols(th.mt.mat.cons.per.mn.2017.gdp ~ unified_rep + DV_lag1 | state + year, data = robust, 
             cluster = ~state)
summary(reg)
```

```{r}
robust <- robust %>%
  mutate(DV_lag1 = dplyr::lag(ff.gdp.perc, n = 1))

reg <- feols(ff.gdp.perc ~ unified_rep + DV_lag1| state + year, data = robust, 
             cluster = ~state)
summary(reg)
```



******************************
***Diminishing Marginal Returns to Climate Policy***


**Do states with higher climate policy have different trajectories?**

*Create a variable for the average level of climate policy over study period*

```{r}
data_climate <- test %>%
  filter(year > 1999, year < 2023) %>%
  group_by(state) %>%
  mutate(climate_policy_mean = mean(climate, na.rm = TRUE)) %>%
  ungroup()

data_climate <- data_climate %>%
  mutate(year_c = year - mean(year, na.rm = TRUE))
```

```{r}
fe_model <- feols(
  Latent_Factor ~ year * climate_policy_mean  |
                  state + year,
  data = data_climate
)

etable(fe_model)

```


**Partisanship**

*Create Partisan Tendency* 

```{r}
# partisan tendency is the unweighted average of partisan share in 
# lower house, upper house, and the Governorship 

# Share of governorship is either 1 or 0 
# Partisan share in legislature is continuous from 0 - 1 

partisan_index <- data %>%
  mutate(
    gov_dem = ifelse(gov_party == "Dem", 1, 0),
    lower_share = two.party.dem.house.share / 100,  
    upper_share = two.party.dem.senate.share / 100
  ) %>%
mutate(
  partisan_index = (lower_share + upper_share + gov_dem) / 3)

avg_partisanship <- partisan_index %>%
  group_by(state) %>%
  summarise(avg_partisanship = mean(partisan_index, na.rm = TRUE))

data <- data %>%
  left_join(avg_partisanship, by = "state")
```


```{r}

fe_model <- feols(
  Latent_Factor ~ year * avg_partisanship  |
                  state + year,
  data = data
)

etable(fe_model)

```


*Do states that begin with high DV values have different trajectories from others? *

```{r}
data_climate <- data_climate %>%
  group_by(state) %>%
  filter(year > 1999, year < 2023) %>%
  mutate(initial_greenindex = first(Latent_Factor[order(year)])) %>%
  ungroup()

model <- feols(Latent_Factor ~ year:initial_greenindex | state + year, data = data_climate)

etable(model)
```


*Is climate policy different in its effect depending on where you start on the DV?*

```{r}
model <- feols(Latent_Factor ~ climate * initial_greenindex + log(gdp.mn.2017.usd) + log(pop.in.th) + manuf.gdp.perc| state + year, data = data)

etable(model)
```




```{r}

df_summary <- data %>%
  filter(year < 2023) %>%
  group_by(state) %>%
  summarize(
    dv_initial = first(Latent_Factor[year == min(year)]),
    dv_final = last(Latent_Factor[year == max(year)]),
    dv_change = dv_final - dv_initial
  )

party <- data %>%
  filter(year == 2020) %>%
  dplyr::select(state, two.party.dem.house.share, avg_partisanship)


df_summary <- left_join(df_summary, party)

ggplot(df_summary, aes(x = dv_initial, y = dv_change)) +
  geom_point(aes(color = two.party.dem.house.share), size = 1.5) +
  geom_smooth(method = "lm", se = TRUE, color = "grey", alpha = 0.3) +
  geom_text_repel(aes(label = state, color = two.party.dem.house.share), size = 3, max.overlaps = 50) +
  scale_color_gradient2(
    low = "red",
    mid = "grey",
    high = "blue",
    midpoint = 50
  ) +
  labs(
    x = "Green Transition Index in 2000",
    y = "Change in Green Transition Index 2000-2022",
    title = ""
  ) + 
  theme_clean() +
  theme(legend.position = "none")

width <- 8

ggsave("Index2000.v.change.index.png", width = width, height = width/1.618, units = "in",
       dpi = 600)


```


```{r}
ggplot(df_summary, aes(x = dv_initial, y = dv_change)) +
  geom_point(aes(color = avg_partisanship), size = 1.5) +
  geom_smooth(method = "lm", se = TRUE, color = "grey", alpha = 0.3) +
  geom_text_repel(aes(label = state, color = avg_partisanship), size = 3, max.overlaps = 50) +
  scale_color_gradient2(
    low = "red",
    mid = "grey",
    high = "blue",
    midpoint = 0.5
  ) +
  labs(
    x = "Green Transition Index in 2000",
    y = "Change in Green Transition Index 2000-2022",
    title = ""
  ) + 
  theme_clean() +
  theme(legend.position = "none")
```



```{r}

climate <- data %>%
  filter(year < 2021) %>%
  group_by(state) %>%
  summarize(
       climate_initial = first(climate[year == min(year)]),
    climate_final = last(climate[year == max(year)]),
    climate_change = climate_final - climate_initial
  )


df_summary <- left_join(df_summary, climate)

df_summary %>%
  ggplot(aes(x = dv_initial, y = climate_change)) +
  geom_point(aes(color = two.party.dem.house.share), size = 1, alpha = 1.5) +
  geom_text_repel(aes(label = state, color = two.party.dem.house.share), size = 3, max.overlaps = 50) +
  geom_smooth(method = "lm", se = TRUE, color = "grey", linewidth = 0.9) +
  scale_color_gradient2(
    low = "red",
    mid = "grey",
    high = "blue",
    midpoint = 50
  ) +
  theme_clean() +
  labs(
    x = "Green Transition Index in 2000",
    y = "Change Climate Policy Score 2000-2022"
  ) +
  theme(
    legend.position = "none",
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 10)
  )


ggsave("Index2000.v.change.cp.png", width = width, height = width/1.618, units = "in",
       dpi = 600)

```


```{r}
df_summary %>%
  ggplot(aes(x = climate_change, y = dv_change)) +
  geom_point(aes(color = two.party.dem.house.share), size = 1, alpha = 1.5) +
  geom_text_repel(aes(label = state, color = two.party.dem.house.share), size = 3, max.overlaps = 50) +
  geom_smooth(method = "lm", se = TRUE, color = "grey", linewidth = 0.9) +
  scale_color_gradient2(
    low = "red",
    mid = "grey",
    high = "blue",
    midpoint = 50
  ) +
  theme_clean() +
  labs(
    x = "Change Climate Policy Score 2000-2022",
    y = "Change Green Transition Index 2000-2022"
  ) +
  theme(
    legend.position = "none",
    axis.title = element_text(size = 10),
    axis.text = element_text(size = 10)
  )



ggsave("Delta.delta.png", width = width, height = width/1.618, units = "in",
       dpi = 600)



```

